Automatically interpreting all faults, unconformities, and horizons from 3D seismic images

Extracting fault, unconformity, and horizon surfaces from a seismic image is useful for interpretation of geologic structures and stratigraphic features. Although others automate the extraction of each type of these surfaces to some extent, it is difficult to automatically interpret a seismic image with all three types of surfaces because they could intersect with each other. For example, horizons can be especially difficult to extract from a seismic image complicated by faults and unconformities because a horizon surface can be dislocated at faults and terminated at unconformities. We have proposed a processing procedure to automatically extract all the faults, unconformities, and horizon surfaces from a 3D seismic image. In our processing, we first extracted fault surfaces, estimated fault slips, and undid the faulting in the seismic image. Then, we extracted unconformities from the unfaulted image with continuous reflectors across faults. Finally, we used the unconformities as constraints for image flattening and horizon extraction. Most of the processing was image processing or array processing and was achieved by efficiently solving partial differential equations. We used a 3D real example with faults and unconformities to demonstrate the entire image processing. Introduction From a 3D seismic image, as shown in Figure 1, one can extract geologic surfaces such as faults, unconformities, and horizons. Faults and horizons are important for seismic structural interpretation because they provide structural maps of the subsurfaces. Horizons and unconformities are important for seismic stratigraphic interpretation because together they construct a chronostratigraphic framework. All these geologic surfaces can be useful for seismic lithology interpretation because they can provide geologically reasonable control for extending a lithology interpretation away from wells. Therefore, extracting these surfaces is a critical part of seismic interpretation. Fault interpretation Automatic interpretation of faults from a seismic image often involves the following four steps: (1) Fault attributes such as semblance (Marfurt et al., 1998), coherency (Marfurt et al., 1999), variance (Van Bemmel and Pepper, 2000; Randen et al., 2001), and fault likelihood (Hale, 2013b) are computed from a seismic image to highlight fault positions. (2) Fault surfaces are extracted from these computed fault images using various methods (Pedersen et al., 2002, 2003; Gibson et al., 2005; Hale, 2013b (3) From extracted fault surfaces, fault slips are estimated by correlating picked seismic horizons (Admasu, 2008) or all seismic reflectors (Hale, 2013b) on opposite sides of fault surfaces. (4) Computed fault positions and fault slips are finally used to undo the faulting in the seismic image (Wei et al., 2005; Wei, 2009; Luo and Hale, 2013). In most methods for fault extraction, the problem of extracting intersecting faults is not well addressed, and the extracted fault surfaces are often represented by triangle or quad meshes, which are often more complex than necessary for subsequent processing. Wu and Hale (2015a) propose to use a simpler linked data structure to represent a fault surface, to extract multiple intersecting fault surfaces without holes at intersection, and to accurately estimate fault slips. Most unfaulting methods (Wei et al., 2005; Wei, 2009; Luo and Hale, 2013) assume that fault geometries need not change to undo the faulting in a seismic image. This assumption makes the unfaulting processing easier, but often results in unnecessary distortions when unfaulting seismic images with multiple faults and, especially, intersecting faults. Wu et al. (2015) propose two methods to simultaneously move fault blocks and faults themselves to undo faulting in a seismic image with minimal distortion. Unconformity interpretation To obtain complete unconformities from a seismic image, we want to extract the angular unconformities with reflector terminations and the corresponding parallel unconformity or correlative conformity with Colorado School of Mines, Golden, Colorado, USA. E-mail: [email protected]; [email protected]. Manuscript received by the Editor 21 September 2015; revised manuscript received 8 January 2016; published online 18 April 2016; Pagination corrected 27 April 2016. This paper appears in Interpretation, Vol. 4, No. 2 (May 2016); p. T227–T237, 14 FIGS. http://dx.doi.org/10.1190/INT-2015-0160.1. © 2016 Society of Exploration Geophysicists and American Association of Petroleum Geologists. All rights reserved. t Technical papers Interpretation / May 2016 T227 D ow nl oa de d 05 /0 5/ 16 to 1 38 .6 7. 12 8. 14 6. R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / conformable reflectors. Most automatic methods can detect only angular unconformities by computing attributes such as seismic coherence (Bahorich and Farmer, 1995) and seismic reflector convergence or divergence (Barnes et al., 2000; Hoek et al., 2010) to highlight areas of reflector terminations. Ringdal (2012) proposes a 2D flow-based method to compute unconformity probability that can highlight termination areas and parallel unconformity or correlative conformity. However, the flow field used in this method is difficult to represent in 3D. For 3D seismic images, this method processes inline and crossline slices separately to compute an unconformity probability volume. Wu and Hale (2015b) propose 3D image processing methods to (1) compute an unconformity likelihood image that highlights the termination areas and the corresponding parallel unconformities or correlative conformities, (2) extract unconformity surfaces from the unconformity likelihood image, and, finally, (3) use the unconformity surfaces as constraints to accurately estimate seismic normal vectors at unconformities and to compute a flattened image with vertical gaps corresponding to the unconformities. Horizon interpretation Automatic seismic image flattening (Stark, 2005; Lomask et al., 2006; Fomel, 2010; Parks, 2010; Wu and Zhong, 2012; Luo and Hale, 2013) is a volume process method to identify all horizons in a seismic image by computing a mapping that transforms the seismic image from the original space into the flattened space. This mapping can be used to extract any number of horizons from the seismic image. However, these methods are unable to match horizons across faults unless additional information such as fault slips (Luo and Hale, 2013) and control points across faults (Wu and Hale, 2015c) are provided. Also, these methods cannot adequately deal with horizons terminated at unconformities unless the unconformity surfaces are provided (Wu and Zhong, 2012; Wu and Hale, 2015b). This paper We present a feasible procedure, as shown in Figure 2, to automatically interpret a complicated seismic image with faults and unconformities, as shown in Figure 1. It is complicated because many intersecting faults are apparent in the seismic image, unconformities, and horizons are dislocated by faults, and many horizons are also terminated at unconformities. To automatically interpret such a seismic image and obtain all faults, unconformities, and horizon surfaces, we apply the following procedure (Figure 2): (1) We first apply the fault-processing methods proposed by Wu and Hale (2015a) to compute fault surfaces and slips from the seismic image and use them to undo the faulting in the seismic image. (2) We then apply the unconformity processing methods proposed by Wu and Hale (2015b) to extract unconformities from the unfaulted image with continuous seismic reflectors across faults and use them as constraints to more accurately estimate seismic normal vectors at unconformities with multioriented seismic reflectors. (3) Next, we use the estimated normal vectors, and use the unconformities as constraints, for a flattening method to compute an unfaulted and flattened image with vertical gaps corresponding to the unconformities. After unfaulting and flattening, one can easily compute horizon surfaces by first extracting horizontal slices in the unfaulted and flattened space and then mapping these slices back to the original space. Figure 1. A 3D real seismic image complicated by faults (red arrows) and unconformities (yellow arrows). Figure 2. Our processing procedure to automatically interpret a complicated seismic image with faults and unconformities as shown in Figure 1. T228 Interpretation / May 2016 D ow nl oa de d 05 /0 5/ 16 to 1 38 .6 7. 12 8. 14 6. R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / Fault processing Because unconformities and horizons are discontinuous at faults, as shown in Figure 1, our first step is fault processing to extract fault surfaces, compute fault slips, and undo the faulting in the seismic image, so that seismic reflectors are more continuous across faults. Fault images and surfaces We assume faults appear as discontinuities that are locally planar in 3D seismic images and linear in 2D slices like those shown in Figures 1 and 3. The key point of this assumption is that when looking for faults, we are not simply looking for discontinuities. Rather, we are looking for discontinuities that are locally planar. With this assumption, we use the method proposed by Hale (2013b) to compute fault likelihood, a measure of locally planar discontinuity, while at the same time estimating fault strike and dip. The fault likelihood image indicates where faults might exist (Figure 3a), whereas the strike and dip images indicate their orientations. In computing these fault images, this method scans over a range of possible combinations of strike and dip to find the one orientation that maximizes the fault likelihood for each image sample (Hale, 2013b). The maximum fault likelihood for each image sample is recorded in the fault likelihood image as shown in Figure 3a, and the strike and dip angles that yield the maximum likelihood are recorded in the fault strike and dip images, which are not shown here. In the fault likelihood image in Figure 3a, we expect relatively high values in areas where faults might exist. However, we do not expect faults to be as thick as the features apparent in the fault likelihood image in Figure 3a. Therefore, we keep only the values on the ridges of the fault likelihoods, and set values elsewhere to be zero, to obtain a thinned fault likelihood image shown in Figure 3b. We also keep strike and dip angles for only the samples with nonzero values in Figure 3b, to obtain the corresponding thinned fault strike and dip images, which are not shown here. These thinned fault images have nonzero values only for samples that might be on faults. As discussed by Wu and Hale (2015a), because most samples in the thinned fault images of fault likelihood, strike, and dip are zero, we can display nonzero samples of the three images, all at once, as fault samples shown in Figure 4a. Each fault sample is displayed as a colored square. The color of each square denotes fault likelihood, whereas the orientation of each square represents fault strike and dip. Fault samples exist only at positions where thinned fault likelihoods are nonzero, and each fault sample corresponds to one and only one image sample. Therefore, these fault samples contain exactly the same information represented in all three thinned fault images. Most fault samples, especially those with high fault likelihoods, are aligned approximate planes, consistent with locally planar fault surfaces. Some misaligned fault samples, often with low fault likelihoods, are also observed in Figure 4a. However, these misaligned samples cannot be linked together to form locally planar fault surfaces of significant sizes. With the fault samples shown in Figure 4a, we apply the method described in Wu and Hale (2015a) to extract fault surfaces (Figure 4b) by linking nearby fault samples with similar fault likelihoods, strikes, and dips. Each fault surface is represented by a set of linked fault samples and not a mesh of triangles or quadrilaterals. The samples on a surface are linked above and below in the fault dip direction and left and right in the strike direction. These links enable us to iterate among seismic image samples adjacent to a fault, in dip and strike directions, as we estimate fault slips. The sets of linked fault samples can be displayed as opaque fault surfaces, as in Figure 4b, by simply increasing the size of each square so that they overlap and appear to form continuous surfaces. However, these surfaces are really just sets of linked fault samples located within the sampling grid of the seismic image; therefore, these surfaces colored by fault likelihoods can also be easily displayed as a 3D fault likelihood image (with nonzero values only at faults) overlaid with the seismic image in Figure 5a. Compared to the thinned fault likelihood image in Figure 3b, spurious fault samples have been removed because they cannot be linked to form fault surfaces with significant sizes. Figure 3. Fault image of likelihoods (a) before and (b) after thinning. The fault likelihoods in panel (a) are displayed in translucent colors. Interpretation / May 2016 T229 D ow nl oa de d 05 /0 5/ 16 to 1 38 .6 7. 12 8. 14 6. R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / We use the fault likelihood image (Figure 5a) to constrain a structure-oriented filter (Fehmers and Höcker, 2003; Hale, 2009) so that it smoothes along structures, but not across faults, to obtain the smoothed seismic image shown in Figure 5b. We then use this smoothed image, instead of the original seismic image, to estimate fault slips. This smoothing does what seismic interpreters do visually when estimating fault slips. We never estimate fault slips using only image samples adjacent to faults; instead, we look at entire fault blocks. This smoothing brings seismic amplitudes from within each fault block up to, but not across, the faults. Fault slips and unfaulting Fault slip is a vector representing displacement of the hanging wall of a fault surface relative to the footwall. In practice, we only estimate dip-slip vectors and relative displacements up or down along the dip of a fault because they are usually more significant than strike-slip displacements in seismic images. Moreover, dip-slip vectors are usually more perpendicular to seismic reflectors, and therefore easier to estimate by correlating seismic reflectors across faults, than are strikeslip displacements. As discussed by Wu and Hale (2015a), to estimate fault dip slip, we first estimate fault throw, the vertical component of dip slip, using the dynamic warping method (Hale, 2013a). Recall that each fault surface in Figure 4b is a set of linked fault samples. The links among fault samples enable us to conveniently iterate on a fault surface and compute differences between seismic amplitudes on opposite sides of the fault. These amplitude differences are computed for every sample on the fault, for a range of shifts that correspond to different fault throws. The fault throws computed by dynamic warping are shifts that minimize these seismic amplitude differences, subject to constraints that the shifts must vary smoothly along the fault surface. Knowing the fault surface with linked fault samples and the computed fault throws, we can walk up or down the fault in the dip direction to compute horizontal components of the dip-slip vector. For each fault sample, the fault throw tells us how far upward or downward along the dip direction we must iterate within the linked fault samples to determine the horizontal inline and crossline components of the dip-slip vector. It is usually difficult to visualize dip-slip vectors in displays of numerous faults. Figure 4. (a) The thinned fault images of fault likelihoods (Figure 3b), strikes (not shown), and dips (not shown) are represented, all at once, as fault samples that are displayed as small squares. Each square is colored by fault likelihood and oriented by fault strike and dip. Links are then built among consistent fault samples, and each set of linked fault samples represents a fault surface that appears opaque in panel (b), where fault samples are displayed as larger overlapping squares. Figure 5. The fault surfaces in Figure 4b are displayed as a fault likelihood image with mostly zeros in panel (a). This fault image is used to constrain a structure-oriented smoothing filter so that it smoothes along structures, but not across faults, to obtain a smoothed seismic image in panel (b). T230 Interpretation / May 2016 D ow nl oa de d 05 /0 5/ 16 to 1 38 .6 7. 12 8. 14 6. R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / Therefore, although estimating all three components of dip-slip vectors, we display in color only the vertical components (or fault throws) on the fault surfaces in Figure 6a. Because the vertical axis of the seismic image is time (not depth), fault throws are here measured in milliseconds, but this method works equivalently for seismic images measured in depth. These surfaces colored by fault throws in Figure 6a can also be displayed as a 3D fault throw image (with nonzero values only at faults) overlaid with the seismic image slices shown in Figure 6b. We observe that fault throws estimated for each fault surface vary smoothly. Also fault throws are non-negative for all fault surfaces, which indicates that the faults shown here are normal faults. As shown in Figure 6, fault slips are estimated on fault surfaces and tell us how to correlate only the samples alongside faults. However, to undo the faulting in a seismic image without distortions, we cannot simply shift just the samples alongside the faults. Instead, we must shift all samples in the image and move entire fault blocks and even faults themselves, as discussed by Wu et al. (2015). Most unfaulting methods (Wei et al., 2005; Luo and Hale, 2013) assume that the geometry of faults does not change and move only fault blocks when unfaulting a seismic image. As discussed by Wu et al. (2015), this assumption usually results in unnecessary distortions when unfaulting a seismic image with multiple faults and, especially, intersecting faults like those on the left of the horizontal slice in Figure 6b. Here, we use the method proposed by Wu et al. (2015) to undo the faulting in the seismic image shown in Figure 6b. In this method, unfaulting shifts skðxÞ are computed in the original space (x 1⁄4 ðx1; x2; x3Þ) by simultaneously solving the following partial differential equations and unfaulting equations: ωðxÞ∇skðxÞ ≈ 0 βcðxaÞðskðxbÞ − skðxaÞÞ ≈ βcðxaÞtkðxaÞ; (1) Figure 6. (a) All three components of fault dip slips are estimated on fault surfaces, but only fault throws and the vertical components of slips are displayed in color on fault surfaces. These colored surfaces can be displayed as (b) a 3D fault throw image (with mostly zeros) overlaid with the seismic image. Figure 7. (a) Unfaulting shifts in vertical, (b) inline, and (c) crossline directions, respectively. Interpretation / May 2016 T231 D ow nl oa de d 05 /0 5/ 16 to 1 38 .6 7. 12 8. 14 6. R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / where ∇ represents the gradient operator and skðxÞ ðk 1⁄4 1;2; 3Þ represent the three components of vector shifts for all samples in an image. Here, ωðxÞ is a weighting function that is zero at image samples adjacent to faults and is one elsewhere. Therefore, this equation is actually defined at all image samples except those adjacent to faults. This equation means that we expect unfaulting shifts to vary slowly and continuously at samples away from faults; tkðxÞ ðk 1⁄4 1; 2; 3Þ represent the three components of fault slip vectors estimated on faults; xa and xb represent samples adjacent to the faults from footwall and hanging wall, respectively. Therefore, this unfaulting equation is defined only for samples adjacent to their faults using the estimated fault positions and fault slip vectors; cðxÞ is a measure of the quality of the estimated slips, and here, we use the fault likelihoods that are already computed on faults; β is a constant number used to balance the partial differential equations and the unfaulting equations. By solving these two different types of equations simultaneously, we obtain unfaulting shifts skðxÞ for all samples in the original space (x 1⁄4 ðx1; x2; x3Þ). These shifts skðxÞ can be easily converted into the unfaulted space ðw 1⁄4 ðw1; w2; w3ÞÞ to obtain skðwÞ as discussed in Wu et al. (2015) because the mapping between the original space and unfaulted space is reversible. Figure 7 shows all three components of unfaulting shifts skðwÞ overlaid with the seismic image. We observe that each component of the shifts is continuous everywhere except at faults. These vector shifts simultaneously move footwalls and hanging walls, and even faults themselves to undo the faulting in the seismic image (Figure 8a) to compute an unfaulted image shown in Figure 8b. In this unfaulted image, all seismic reflectors are well aligned across faults and no distortion is observed. Unconformity processing From an unfaulted image (Figure 8b) with more continuous reflectors across faults, interpretation of seismic horizons and unconformities becomes more straightforward. Before extracting seismic horizons, we want to first interpret unconformities because horizons might be terminated at unconformities and the interpreted unconformities can serve as constraints for manual and automatic horizon extractions. To obtain complete unconformities from a seismic image, we want to extract angular unconformities with reflector terminations and the corresponding parallel unconformity or correlative conformity with conformable reflectors. Most automatic methods (Bahorich and Farmer, 1995; Barnes et al., 2000; Smythe et al., 2004; Hoek et al., 2010) can detect only angular unconformities. Here, we use Figure 8. A seismic image (a) before and (b) after unfaulting. Seismic reflectors are more continuous across faults after unfaulting as highlighted by the red arrows. Figure 9. Two seismic normal vector fields (a) us;c (green) and us;a (yellow) are computed from structure tensors that are constructed with vertically causal and anticausal smoothing filters, respectively. These two vector fields are different near the unconformity, and (b) unconformity likelihood is an attribute that evaluate the differences of the two vector fields. T232 Interpretation / May 2016 D ow nl oa de d 05 /0 5/ 16 to 1 38 .6 7. 12 8. 14 6. R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / the method proposed by Wu and Hale (2015b) to compute an unconformity likelihood image that highlights the termination areas and the corresponding parallel unconformities, as shown in Figures 9b and 10a. As discussed by Wu and Hale (2015b), the unconformity likelihood is defined as an attribute that evaluates the differences between two seismic normal vector fields estimated from a same seismic image. These two vector fields correspond to the largest eigenvectors of the two structure-tensor fields Ts;c and Ts;a constructed from the same seismic image but with different smoothing filters: Ts;c 1⁄4 hhg1g1isic hhg1g2isic hhg1g2isic hhg2g2isic (2)